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ABSTRACT 

An Implicit finite difference procedure is developed to solve the 
unsteady full potential equation in conservation-law-form. Computational 
efficiency is maintained by use of approximate factorization techniques. 

The numerical algorithm is first order in time and second order in space. 

A circulation model and difference tions are developed for lifting 
airfc s in unsteady flow, however, thin airfoil body-boundary conditions 
have been used with stretching functions to simplify the development of 
the numerical algorithm. 

1. INTRODUCTION 

Current numerical algorithms to compute unsteady transonic inviscid 

flow about complex configurations are frequently either inadequate or 

too costly to use for routine analysis of a large class of two and three 

dimensional flowfields. In particular, numerical algorithms* for 

unsteady transonic small disturbance theory neglect terms that can be 

important, for example, in helicopter rotor flowfield simulation. Numerical I 

I 

algorithms based on the Euler equations are suitable for any inviscid [ 

flowfield simulation, but current numerical algorithms for the Euler 
equations have large computer time and computer storage requirements. The 
successful development of -^n unsteady conservation-law-form full potential 
finite difference algorithm therefore offers a practical design tool. 

Unsteady potential theory can satisfactorily replace the Euler equation 
solutions if the shock waves are sufficiently weak and if the equations 
employ an equivalent circulation model. Yet, the full potential equation 
has similar computer requirements in time and storage to the simplified 
small disturbance theory. 
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The purpose of this paper Is to present an approximate-factorization. 
Implicit finite difference algorithm for the unsteady full potential 
equations In conservation-law-form. Conservation-law-form Is maintained 
In order to "capture" shocks with correct strength and location. During 
the course of this research, two other unsteady procedures have been 
developed for the full potential equations in conservative form. One 
procedure, due to Chipman and Jameson" solves a system of three equations 
for three unknowns. The other, due to Goorjian,* is a scheme which is 
quite similar to the one developed independently here. In this report the 
governing equations are first reviewed in Section 2. Differencing of the 
governing equations in conservation-law-form and their numerical solution 
are discussed in Sections 3 and 4. Finally results and concluding remarks 
follow. Throughout small disturbance boundary conditions are used to 
avoid unduly complicating the development of the numerical method for the 
full potential equations. 

2 . GOVERNING EQUATIONS 


Stretching T ransfo rmations 

As governing equations for two-dimensional unsteady full potential 


flow we use the conservation of mass equation 


^ ^ ^ 
3t 3x 3y 


= 0 


( 1 ) 


where density, p, is determined from the unsteady Bernoulli relation 

1 2 

y 


D = {1 + (Mf, - 2^^ - 


V ^ / 


Here = u, 


V satisfy the condition of irrotational ity 


3u ^ 
3y * 3.x 


0 


(3) 
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and in deriving Eq. (2) the far field flow is assumed to be steady. 
Throughout the equations have been referenced to freestream quantities, 

p » p/p^ , u » u/a^ , X « x/il 

and the tilda has been suppressed. Here a is the sound speed while 
£ is a reference length such as airfoil chord, c. 

Boundary conditions are best satisfied by mappings which can also 
be used to cluster grid points so as to en-tance numerical solution 
accuracy. As shown by Viviand® conservation-law-form of tq. (1) can 
be maintained under the general transfornation 

i = f;(x, y. t) 

n = n(x, y. t) (4) 

\ = t 


gi ving 

^-(P/J)+ (pU/J) + 1^ (pV/J) = 0 

where J is the Jacobian ■' - i r\ and 

'X y y X 

u = f-t * - i) * - y) 

V ■ Ht * - i) * riy(4'j, - y) 

The terms 4> and are expanded by chain rule 
X y 

C' = f + n 't’ 

^x ^x^f. X n 


(5) 


( 6 ) 


(7) 
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and U and V can also be written 


U = ^ + 4>^Vn) = VC • (^^vc + (t.^Vn - r^) (8a) 

V = + Vn • (ij^vc + 4i|^Vn) = Vn ‘ + '{'^Vn - f^) (8b) 

-► f 

where r = (x, y) . Subject to the same transformation the Bernoulli 


equation is written 


t«i - 2(*^ . - (5^ * C^)*^ 


2(C n + C n )<t>r<i> - (n^ + 


For the present purpose of demonstrating an unsteady full 
potential algorithm in conservation-law-form the transformation is 
restricted to the form 


^ = C(x) 

( 10 ) 

n = n(y) 

and the thin-airfoil body-boundary approximation is utilized. Because 
of the transformation Eq. (lO), the governing equations then simplify to 


((p/JK^«j) * I;; ((P/J)n^) = 0 


P = (1 + 


- (M„ - - HyO^)} 


Boundary Conditions 


As boundary conditions, tangoncy (i.e. V = 0) is imposed on 
the body surface and tie flow at infinity is required to be uniform and 
steady. The thin-airfoil body-boundary approximation is utilized 
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since our sole interest here is that of developing an efficient unsteady 

algorithm. The complete mapping of boundary conditions unduly complicates 

this task. Surface tangency (il' - x)n + (({> - y)n * 0 is replaced by 

X X y y 




y 


u ^ 

body slope 


( 12 ) 


which is imposed on the grid lines adjacent to the airfoil -- see 
Fig. 1. 

At distances far from the body u = u and with our use of a as 
a reference velocity 



where M is the freestream Mach number. Here we chose to cant the 

OP 

airfoil to obtain angle of attack rather than have the flow vector be 
specified at an angle as is often done. 

‘■'^r lifting cases we must allow for a jump of potential across 
a cut. As we intend to partially model an unsteady shear layer leaving 
the airfoil, this cut is taken off from the trailing edge, see Fig. 1. 

On the back boundary then. Eg. (13) does not apply and so on the bacr. 
boundary freestream pressure (or density in isentropic flow) is imposed. 
From Eq. (11b) this implies that 

^ ( 14 ) 

along the back face cd of Fig. la. Neglecting small perturbations, this 
simplifies to 

(t + f. M (f = (15) 

^ T X. <*• L ^ 

Eq. (13) is imposed along ail other far field boundaries. Fig. la. 
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In unsteady flow, vorticlty Is continuously shed from the 


In a potential flow formulation this phenomenon can be modeled with a 
cut aligned with the shear layer. Here we further Idealize this cut by 
keeping it In the mean chord line plane as sketched In Fig. lb. 

The jump In potential along the cut, A<t = r, will vary with time, 
by Imposing the usual shear layer assumptions an equation for r(x, t) 
along the cut can be derived. Across an Inviscid shear layer normal 
velocity and pressure (p) are continuous, while tangential velocity 
and density can be discontinuous. Here, however, the incoming flow is 
uniform Isentropic (i.e. homentropic) and the shocks are weak. Thus, 
we have previously assumed that p = throughout, and density Is also 
taken to be continuous across the cut. 

Define r = (I 41D " 4 ’^ " across the cut as illustrated in 
Fig. 1b. From the condition that the velocity normal to the shear layer 
Is continuous. Hilly]) = 0 

Vt " Vu " ( 16 ) 

The Bernoulli relation together with the continuity of density and normal 
velocity require., that 

(20^ + 0^ 0y)^ = (20^ 

or 

( 2 »J * ♦ r) + [3^(»u ♦ r)]^ 

* * r)]^ 

SO that 
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Consistent with thin airfoil body boundary conditions, Eq. (17) Is 
approximated by 

Tt + MJx = 0 (18) 

Eq. (18) gives an equation to determine r(x, t) along the cut. 

Although is continuous across the cut and * 0, in unsteady 
flow and 4 )^^ is discontinuous. This follows immediately from 

which on employing density continuity becomes 

C3^(P0^) + P3y 0yll = 0 (20) 

or alternately 

Ptt-fyy]] = - 3j<(pr^) (21) 

This last equation will be needed to develop accurate difference 

expressions for Eq. (11a) when it is differenced to just either side of 

the cut (see Section 4). For steady flow = 0 and from Eq. (17) or (18) 

it follows that r = 0. Thus, in steady flow <t> is continuous across 
X yy 

the cut. 

3. CONSERVATIVE DIFFERENCING AND LOCAL LINEARIZATIONS 
Ti me Derivat iv e Term 

Stable conservative difference operators have been successfully 
developed for the steady state part of Eqs. (1) or (11a) (see Refs. 7-10) 
so the main task initially undertaken in this research effort has been 
to develop stable conservative difference operators for the unsteady 
terms as well. Development of nonconservative time differencing schemes 


7 


has always appeared to be straightforward. Because p * p(4s), Eq. (1) 
can be rewritten as 


3o 94> ^ 
3$ ot 



♦ ^ 


9p M 
X 30 3x 


3d> By 


where 


3 £ 

3^ 





0 


( 22 ) 


(23) 


is a differential operator that does not cormute. Performing the indicated 
operations, Eq. (22) is found to be a second degree wave equation in 


nonconservative form 


♦tt * ■ (p"'-' - * (p'-' - (24) 

1 2 

here using p^" = a^ (in the chosen nondimensional variables) gives thn 
familiar textbook form of the potential equation (see, e.g. Ref. 11). 
Various finite difference schemes have been advanced for simpler but 
similar forms of this equation (1-3) and we assume they can be made to 
work for Eq. (24). 

The problem arises, though, that Eq. (24) cannot be brought back 
into conservation-law-form with (J) retained as the dependent variable. 

If, however, one is willing to introduce into the differential equations 
the same expansions used in deriving the difference formula, conservative 
form can again be maintained. By Taylor expansion 


' “o * 


(25) 


where o represents a neighboring known state or solution. Substitution 
of Eq. (25) into Eq. (1) gives 
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(26) 



(P4>y) 


Assuming that all of the density terms are formed with exactly 
the same difference operators, we assume that the local linearizatio.. term 
^t (*^o ‘ If^o '^o) ^ higher order contribution to numerical 

stability. Moreover, if 0 - i||q is small, for example or 

()>. - (ji. , where t * nAt and x = jAx, the error due to expanding p is 

J J * 

second order accurate and is no greater than that usually made in time 
di fferencing. 

Consequently, Eo. (26) is a conservation-law-form of Eq. (1) 
which is assumed to have linear stability properties equivalent to the 
equation 


3 t ^“ 1 *^ 


t " “ 2 ^ “3^y^ ' 


(27) 



I 

'i 


As with the nonconservative form, Eq. (24), one can devise stable 
difference schemes for Eq. (27) when , a^* <^ 3 * P considered 
as constants. With care, these schemes can be expected to apply for the 
nonlinear equations as well, although, as discussed immediately below, 
assuming p to be a constant fails for transonic flow. A correct treatment 
of the right-hand side of Eq. (27), however, is known from Refs. (7-10). 


i]|. 

1 
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space Deri .a lives 

Just as the local linearization approach guides the proper 
treatr^nt of the time derivative term, it can be used to analyze the 
space differencing schemes developed for the spatial terms of Eq. (1) by 
Jameson,’ Holst and Ballhaus,* and Hafez, South, Murman‘® for transonic 
flow regions. Before proceeding with this analysis, however, some 
prelimimry background information is in order. 

Consider first the simple steady state model problem 

0 - * ‘yy ' » 

and let V , a , V and A be defined in the conventional way, e.g. 

X X y y 

' (*k+l ’ known that the 

following difference schemes 


(1 - Vy*^ " ° 

M < 1 

OD 

(29a) 

( 1 - + VyAy0 = 0 

M > 1 

OD 

(29b) 


are respectively suitable for the model elliptic (M^ < 1) and the model 
hyperbolic (M^ > 1) problems. The differencing (29b) is divergent if 
< 1 while that of (29a) is convergent for > 1 only if 
|Ax/((l - wf^)Ay)l < 1, an impractical constraint for -* 1 . Murman and 
Cole,*’ aware of these simple constraints, introduced their type 
dependent different operators for the transonic small disturbance 
equation. In this approach the streamwise spatial difference operator is 
switched from central to backward as the local Mach number is less or 
j, '"eater than 1 . 
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Various refinements to the Murman-Cole differencing have bt.n 
introduced, but apparently overlooked is the fact that Eq. (28) can also 
be differenced as 

• ”^x'’x» * 'yV ■ “ 

without any need to switch difference operators. An extra boundary 

condition is needed in the x-direction (which is also true of many higher 

order difference schemes) and one must pay attention to the physics 

2 2 

because replacing by is also convergent for any The 

differencing scheme given by Eq. (30) is also less accurate than Eq. (29a), 
but, subject to the usual approximations such as use of periodic boundary 
conditions, the authors can demonstrate the convergence of Eq. (30). 

Although apparently not recognized as i>dc.i, the di.'ference 
schemes of Refs. (7-10), which are used when the local Mach number is > 1, 
actually mimic the difte> crcing scheme represented by Eq. (30). The 
same differencing is not used in subsonic regions because the pure 
central differencing scheme is more accurate. The local linearization 
process will now be used to clarify these remarks. 

Consider the spatial derivative term Use of the local 

linearization, Eq. (23), with this term gives 

»x»«x ■ * (»xlf * 

■ >xt‘=o'''x» ■ <»x'>^‘’')o<’t * ‘x'o’x ■* »yl<,5,X« - «0>1 
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The and terms of Eq. (31) are not crucial to our discussion, 
and we can justify ignoring them by assuming steady small perturbation 
flow. The bracket term of Eq. (31) is clearly a type dependent term of 
the form 

3 3 (b - 3 p (u /a )^3 (t (32) 

x^o x^o' 0 o' x^ ' ' 

where the subtilda is used as a reminder tha'. this operator is the tor 
used in the Bernoulli equation to form p. U;'.ing the differencing formula 
of Refs. (7-10), the term of Eq. (1) i . always treated by central 

difference formula, although in supersonic flow regions p is evaluated 
shifted back to the previous point in x. That is, the term 
replaced by (see Ref. 8 or 9) 

f V = 0 ,(u/a)^ < 1 
( ? ( 33 ) 

^ ^ ^ ^ = l,(u/a)^ > 1 

where p is a shifted backward (i.e. upwind) value of density. The 

precise difference operators are treated in Appendix A; but, the crucial 

point to be made here is that when density ’s shifted backwards, the 

operator ? of Eq. (32) also is shifted backwards (upwind) as this 
/\ 

/^✓ 

operator comes from the local linearization of the Bernoulli equation. 

Thus the differencing Eq. (33) is equivalent to the differencing 

6jj[(l - v)p + vp]6j^<() = 'SxPo^x'^ ■ ‘^x^o(“o^®o^'^x‘^ + . . . (34) 

central central v = 0 
backwards v = 1 

where 6^ are symbolic difference operators. A rigorous discussion of 
the actual difference operators is given in Appendix A. 
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If the local linearization argument given above is indeed valid, 
then the differencing represented by Eq. (34) with v = 1 should be 
stable in subsonic as well as supersonic region so that switching is 
unnecessary. This is indeed the case as reported in Ref. 8, and, with 
the use of second order differencing, the unswitched differencing scheme 
gives quite satisfactory results as we later demonstrate. 

4. NUMERICAL ALGORITHM 

An implicit approximately factored (AF) finite difference scheme is 
developed for Eq. (11) in conservation-law-fori-» by using the local 
linearization discussed in the previous section. The advantages of an 
implicit AF algorithm have been discussed before, see, for example 
Refs. 1, 2, 8, and 13. The algorithm presented here is optionally first 
or second order accurate in space and time, but the second order in 
time option has not proved to be satisfactory. 

An implicit finite difference scheme for Eq. (11a) is given by 
(here suppressing spatial indices for convenience) 


/.,R+1 ..n, ^ r/r2^\n+lT- .n+K . -t r, 2^,n+lx aI+It 

(p - p ) + h6^[(C^p) 6^(J) J + h5^[(nyP) 6^(t) ] 


= (a/3)(p" - a'’’^) 


(35) 


where p i (o/J), h = a = 0 for first order time accuracy 
or a = 1 for second order time accu*"acy. The operators 6^ and are 
defined at each j, k point (x = jAx, y = kAy) by 
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(1 + 0)p. + (1 - 0)p. , 1 

* ^ - ♦j) 

p [ Pi + Pi 1 (1 + 0)p- 1 + (1 - e)p. 

■ ^V'^^j-1/2^^ ■ 2"^+ Vj J ^ ■ '^^j-lH36d) 


■ ^^y^‘^^k-1/2 2 ^ ^*^k ’ *^k-1^ 

here AC = An = 1 and only the varying indices are indicated. The 
parameter 6=2 for second order spatial accuracy in supersonic regions, 
and 0 = 1 for first order accuracy. The switching parameter v is defined 
in a way similar to Refs. 7-10 and 

V = [1 - (p/p*)^]c 1 < c < 10 (37) 

V E 0 if V < 0, i.e. subsonic 

V E 1 if V > 1, i.e. supersonic 

The parameter v can be set to one throughout, but accuracy will be 
impaired unless e is also set to 2. The operators (36a) and (36b) 
assume that the flow will be supersonic only in the positive x-direction 
(see e.g. Ref. 9 for extensions). The density is found from the Bernoulli 
equation with (AC = An = 1) 


(38a) 


" Vj 

^Ik = " Vk 

<l>^r'’^ = [(2 + a)(<^"^^ - (t") - a(^" - <f"‘b]/(2At) = (38c) 

The metrics £ and n are obtained from 

y 

«x = 2/(Xj„ - Xj.,) (39a) 

^ = - >'k.l> <35l>) 

2 

while the term j+1/2 formed either as 

(5x/J)j.l/2 ” " <«x/J)j]« (^0^> 

or 

2 2 2 

The terms j-1/2’ ‘'^y'^‘^^k+1/2’ ^’^y'^‘^^k-1/2 similar treatment. 

If Eq. (40a) is used, it is essential to add 

to the right-hand side of Eq. (35) to subtract out a numerical truncation 
error due to incomplete metric cancellation (see also Refs. 14 and 15). 

The term, relation (41), should equal 0 where = M^x but this is not 
obtained if Eq. (40a) is used and the error can be very appreciable on a 
highly stretched grid. A similar correction in the Bernoulli equation 
isn't needed due to the choice of difference operators and stretchings. 

The local linearization of Section 3 are now introduced. Using 
Eq. (25) to both expand about and to expand p^ about , the 
terms p^^^ - p^ of Eq. (35) are approximated 

1 5 


- e" . (p-'/j"*’ - (8"/j"*')[6, . (£2)"*' <6^](r' - ♦'')) 

- /p"-'/j" - (8"‘'/j")[s, t - ♦"■’)') 


where of Eq. (25) is replaced by Eq. (23) and 8 = p^'^. Here also the 
coefficients (f- and 4> are evaluated as indicated by Eq. (38), while 6 , 6j., and 6 

s n X ^ 

are the same operators defined by Eq. (38). Note both and p^ are linearized 
not just p" , otherwise, conservation-law-form and second order accuracy is lost. 
The operators 5^, 6^, and 5^ were defined by Eq. (38) because their form is 
dictated by the selection of difference operators used in the Bernoulli equation 
for density. These operators work on any product of elements to their right. 

The second derivative terms are rear.anged into delta form and p*^^^ 
is eliminated as follows: 


T / \ ^ T" 


6,(f2/J)"^\.%(r' - 0 ") 

+ + a(p" - p"‘^)]^^4>" 


(43) 


The n-term is treated in a similar fashion, and for a = 1 second order time 
accuracy is maintained. 

Applying Eqs. (42) and (43) to Eq. (35) gives the locally linearized 
form of the implicitly differenced governing equation: 
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- hZ^(n^/0)"*'p\)(r' - ♦") 

= (b"‘’/j")[ 5^ + («^"'»'r'*5 ^ - 6"'') 

t h(j5(t^/J)"*’[p" * p(p" - + T^(n^/J)"*'(p" + a(p" - p"'')F^j (44) 

where one must remember that the operators work on any product of terms to 
their right, e.g. 

(6^a6^ + 6^66^) (c + d) = <S^a6^c + *5^a6^d + 5^b5^c + <5^b6^d 

The operator 6^ that appears in Eq. (44) is given by Eq. (38c). As 
it operates on A(^, however, we chose to replace it by the operator 

5/^' = - <l'")/(At) = (1/At)(I - E■’)<^"^^ (45) 

which agrees with Eq. (38c) only if a = 0. Eq. (45) is now substituted into 
Eq. (44), e'^/( At) = (o^'^)'^/( j'^^^At) is divided through, and the left-hand 
side is approximately factored into £; and n operators: 

(I + At(nJ)"^’0"6^ - At(/^V6'’)h7^(nJ/J)"*V'V^}x 

{I + At(E^)""’<j>"6, - At(j"^V8")hI^(C^/J)"^^p'^^}(^"^^ - <f") 

= [1 + (6"'Vb")(J'’^Vj")](4" - - (b"'Vb")(j"''Vj")(4>"’’ - ^"'^) 

+ At(B"’'VB")(/"Vj")[(C^)^J'’6^ + (ny)%";^6^](4" - 

+ At(/'VB"){[1 - (cr/3)](b" - p""b + hJ^(C^/J)"*^[p" - a(p" - 

+ hX^(ny/J)'^^^[p^ + a(p'^ - p^’^ )]7^4>'^} (46) 


17 


This equation has the form 




and it is implemented as an algorithm as 


= R 


(48a) 


L^A(|> = A(}>* 

X Aa" 

(P = (p + a<p 


(48b) 


(48c) 


The algorithm Eq. (48) requires only a series of scalar tridiagonal 
inversions and it is therefore very efficiently implemented. Computer 
storage equivalent to four levels of have to be supplied with p computed 
from the Bernoulli equation as needed. Alternately three levels of (Ji and 
one level of p can be used, although in our test program two levels of p 
where stored for convenience in programing ' computer code that changed from 
day to day. All exponential functions were eliminated by using binomial 
expansions of the form: 

P = 1 + <pe(] ^ eO + ^ e))) (49) 


where 


E = (y - IXK^ - 2*, - - nyt^)/2 

♦ ' 1/(y - I) 


6 • s*(l E *£(1 ♦ c(l * . . .))) 


J 


ty 

where the constant 6* " is evaluated using the nondimensional 

critical density, and 

c * (p/p*) - 1 
fp • 2 - y 

The expansion for density gives at least four-place accuracy in the flow 
regime of interest. The expansion for 3 need not be so accurate as it is a 
linearization coefficient. On the CDC 7600 a run time of 90 seconds was 
reduced to 60 seconds because of the use of the expansions. In this example 
the solution was advanced 400 times steps on a 50 x 62 grid. 

Implementation of Boundary Conditions 

The thin-airfoil body-boundary tangency condition, Eq. (12), was 
directly built into the implicit algorithm. Indeed, it was in order to 
satisfy this boundary condition that n-tridiagonals are formed before the 
E-tridiagonals in Eq. (46) or (47). The difference equation (+/- indicate 
upper or lower surface) 



is implemented into the numerical algorithm as follows: Eq. (48a) is 

first formed just as before ignoring the boundary. Then for points 
immediately above and below the airfoil chord -- which is centered in 
between, see Fig, 1 -- Eq. (51) is used to correctly overload the tridiagonal 
elements as sketched ir Fig. 2. Likewise, in forming the E-tridiagonals 
adjacent to the body-boundary, elements are first loaded as usual and then 
those elements corresponding to points just above or below the body are 
overloaded as sketched in Fig. 3. This approach insures that A(p above and 
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below the airfoil are unaltered by the Xl-inversions, yet come in Implicitly 
for Xl-derivati ves just ahead and behind the airfoil leading edge and 
trailing edge. 

Adjustments to the difference equations and numerical algorithm 
must also be made for points just above and below the cut indicated in Fig. 1. 
For these points, is altered as follows (see Appendix B) 

^Ik = (Vl ■ "^k-l ■ ^ [J/(16pny)]I^(pC^/J)y (52) 

where this change effects the calculation of density as well as coefficients 

in Eq. (46). As indicated by Eq. (21), the difference formula for 
2 

must also be adjusted for points adjacent to the cut. For points 
above the cut (see the derivation in Appendix B) 

<V''>ktl/2 2— <Vl ■ »k> 

- (♦,, - ♦u., - r) - 1/8 ?j(pC^/J)T^r (53) 
while for points below the cut 

' <''/''>k.l/2 '♦k*! - *k - '■) 

- (ly/J)|,.|/2 - Vl> * 1/8 'S’j(pE^/J)rjr (64) 

2 

For simplicity the term 1/8 is lagged at time level n and it thus 

only enters into the right-hand side of Eq. (46). 

At the end of each update of the field for (pt new values of r are 
obtained along the cut by solving Eq. (16). This equation is differenced as 




.n+1 


- ‘S-lx’'"*' 
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with Initial data in x supplied from the known value of at the 

airfoil trailing edge just upstream of the cut. 

5. RESULTS 

Calculations to verify the accuracy of the numerical algorithm are 
shown in Figs. 4 to 7. The Cp distribution for steady-state flow about a 
6 percent thick biconvex airfoil is shown in Fig. 4 for M = 0.857 and 
a = 1“ angle of attack. The conservative full potential result obtained 
using the second order option, 0 » 2, is compared to experiment and a small 
disturbance nonconservative result. The conservatively captured shock is, 
of course, downstream of the nonconservative result. In this case the 
vertical far field boundaries are placed at 20 chord lengths away from the 
body, although equivalent results are obtained if these boundaries are as 
close as 12 chords. 

As shown by the results of Fig. 5, essentially the same steady-state 
solution is obtained if the switching parameter v is always set equal to 1 
(i.e., no switching), provided the second order differencing 9 = 2 is 
used. Slight differences in the two solutions are observed at the leading 
edge singularity and at the shock wave. The unswitched scheme, v = 1, 0 = 2. 
gives a more pleasing shock wave result. 

Unsteady solution accuracy is demonstrated by a comparison with a 
small disturbance low frequency result obtained by Ballhaus and Steger.* 

In this test case, the airfoil is a biconvex profile which varies in time 
from zero thickness to 10 percent thick, and then thins back to zero. A 
shock initially forms past mid-chord, and then moves back to the trailing 
edge as the airfoil thickens. As the airfoil then thins, the backward 
motion of the shock stops, and then the shock propagates upstream. Fig. 6 
shows a trace of mid-chord pressure for such a case with a comparison to 
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the low frequency small disturbance result. The small discrepancy between 
the theories is about what one would expect for this case, and similar 
results are obtained with the full potential schemes of Goorjian.* 

As a final test case the unsteady flow about an impulsively plunged 
flat plate was computed with * 0.8 and a plunge velocity equivalent to 
0=1°. Linear theory is valid for this case and good agreement with linear 
theory is obtained by Chipman and Jameson” who solved the full potential, 
as well as by Beam and Warming*® who solve the Euler equations. Good 
agreement with linear theory is also obtained with the present code as 
the load distributions plotted in Fig. 7 illustrate. 

6. CONCLUDING REMARKS 

An implicit finite difference procedure was developed to solve the 
unsteady full potential equation in conservation-law-form. Local 
linearizations were successfully used to derive a correct time differencing, 
avoid iterative between time levels, and to correctly analyze the spatial 
differencing. An unsteady circulation model was developed and various 
test cases were done to verify the accuracy of the numericv’l algorithm. 
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APPENDIX A 


Details are given In this appendix to show, to first approximation, 
the equivalence of Eq. (31) for the actual difference operators used In 
the Holst-Ballhaus delayed density scheme. Holst and Ballhaus have 
previously shown the equivalence of their scheme to the artificial 
viscosity method of Jameson’ and of Hafez, South, and Murman. ” 

To make the steps clear, the analyslb will be detailed here for the 
fully delayed case, v = 1 In Eq. (33), In one-dimensional steady flow. 

Then to first approximation It will be shown that 



where 
e.g. = 


In the Bernoulli equation Is computed with the midpoint rule. 



(A. 2) 


To verify Eq. (A.l) the terms are expanded as follows: 


('J'j+1 ■ 

‘^j-1/2 "^Ax - Pj-l/2\‘^j+l/2 


1 


* [p (Pj-i/2 ■ ♦ '^j+l/2 ■ ^ 

+ (Pj-1/2 ■ 

“ ^^x*^j+l/2 ^x^^j-1/2 ■ 

where ? and p are from a known nearby solution and 
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From the steady one-dimensional Bernoulli equation 

1 

e • [1 ♦ 


(A. 4) 


and via Taylor series 


p » p - - 'I') O(A^) 


(A. 5) 


Using (A. 5) to eliminate Pj_i/2 ‘ ^ gives 


‘*J4’ ' *j’ . 


J-I/2 - ^ ‘ oVjtt/2 - ■ ♦> ’ 


(A.C) 


where the midpoint rule Eq. (m 2) is again used to approximate derivatives. 
Likewise 


"j-a/J ■ AVj.,;2 - 


(A. 7) 


Subtracting (6) from (7) and dividing by Ax gives Eq. (A.l) plus a term 
of [O(A^) - 0(A^)]/Ax - O(A^). 

Other ways of shifting the density function give the first 
approximations : 

(centered) 

- * j! • (A. 8) 

Ax^ * AX^ 


(centered - upwind) 

(2Pj., 22 - “j-3/2’l*Jtl ■ ‘j’ ■ ’^‘'j-3/2 ■ '’j-6/2"*j ' *J-l’ 


. i ' '-y * -8-3-2 »'».1 • - “^.i-2 * ♦j-3’ 

both of which are second order accurate approximations. 


(A. 9) 
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The difference schemes used In Section 4 do not use the mid-point 
rule to evaluate In the Bernoulli equation. Rather Pj_]/2 
approximated by (Pj + pj .|)/2 and Is approximated with the three point 
relation, Eq. (38a). This gives the first approximation 


(p. + P.-_-|) (p. I + P.p) 

— J — .-II- /A. . - A.l J 


T 


- *. 1 ) - 




Ax 


= - . p2-Y^2 ("PjAi - »j-3^ 


Ax 


(2Ax)‘ 


(A. 


Unlike Eq. (A.l) a pure upwind formula is not used with the negative 
2 2 

coefficient, . The resulting formula, however. Is sufficiently 

biased to obtain the same effect. In fact, the combined differencing 
can be more upwind than (A.l) since the coefficient to the j+1 point Is 
p(l - u^/(4a^)) rather than 'p. 


APPENDIX B 


The derivation of Eq. (53) Is sketched below using the notation 
indicated in Fig. B.1. The derivative 3 p3 <|> is first approximated as 


^ ^ my? 


o(Ay) 


Now 


\ ^ ^ *^y|u * (t) <^yy|u 




k-1 




*yy|£ 


or 


^y|k-l/2 

'^ylk-1/2 


•JiyL = (Jiy u is found to be 


= - Vi) - r 


(Ay)‘ 




~E~ u 


Ik')' 


<^yy|o))/Ay 


or 


^^'^y|k-l/2 “ ^“^k ■ “^k-l ' ^*^yy^ 


App>*oximating Eq. (21) as 


^k ^ ^\-l 


Pk " Pk-1 


2 ^ yy“ X 2 x 


and substituting Eq. (B.2) and Eq. (B.3) into Eq. (B.l) gives 
. ^^k+1 ■ ^k^^'^k+1 ■ *^k^ ■ ^^k ^ ^k-l^^'^k " \-l " 


3yP3y>^ ^ 


2{Ay)‘ 


. 1 T /'k ^ ^'k-l’lT . 

" 5 2 


(B.l) 


(B.2) 


(B.3) 


(B.4) 
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Eq. (53) differs from (B.4) only in the addition of metrics. 

The derivation of Eqs. (52) and (54) proceed in a similar way where 
in deriving Eq. (52) one starts with (for k above the cut) 
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(a) Tridiagonal Lf^A4>* = R before tangency 
boundary conoition 



(b) Tridiagonal = R with body tangency 


Note: Indices on A, B, C, and * on A4>* are deleted; 

' hower- " ' '‘upper 


Fig. 2 . Adjustment of n-tridiagonal to implicitly include 
the thin-airfoil body-boundary condition 
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(a) Tridiagona! before adjustment 

for tangency boundary condition. 


ABC 

ABC 

1 

1 


= A*}** 


ABC 


(b) Tridiagonal adjustment if and only if 
^ * *^lower *^lower 


Fig. 3. Adjustment of C-tridiagonal along k 
■ "upper- 


lower 
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CIRCULAR ARC PROFILE, r » OM 
M„- 0.857, a-1“ 

A ▼ EXPERIMENT. NASA TN D-16 

-.8r NONCONSERVATIVE. SMALL DISTURBANCE 

UNSTEADY FULL POTENTIAL 



0 1.0 

x/c 

Fig. 4. Steady state flow about biconvex airfoil 
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x/c 


Fig. 5. Comparison of solutions obtained with 0 « 1.8 in Eq. 36a 










(b) Solution at 0.8 chords of travel 


Fig. 7 . Continued 



FLAT PLATE 


LINEARSOLUTION 

POTENTIAL SOLUTION 

BEAM WARMING EULER SOLUTION 


(c) Solution at 2.4 chords of travel 



q. 7 . Continued 



